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Abstract 

Proper orthogonal decomposition methods for model reduction utilize information about the solution 
at certain time and parameter points to generate a reduced space basis. In this paper, we compare 
two proper orthogonal decomposition methods for reducing large systems of ODEs. The first method is 
based on collecting snapshots from the solutions only; the second method uses snapshots from both the 
solutions and their time derivatives. To compare the methods, we derive new bounds for the 2-norm of 
the approximation error induced by the each of the methods. The bounds are represented as a sum of two 
terms: the first depends on the size of the first neglected singular value while the second depends only on 
the spacings between the snapshots. We performed numerical experiments to compare the errors from the 
two model reduction methods applied to the semidiscretized FitzHugh-Nagumo system and investigated 
the relation between the behavior of the numerically observed error and the error bounds. We find that 
the error bounds, though not tight, provide insights and justification for using time derivative snapshots 
in POD model reduction for dynamical systems. 
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1 Introduction 


In many areas of science and technology, complex multi-physics time-dependent problems are modeled by 
large systems of differential equations. Their analysis often poses huge computational challenges as it requires 
multiple, prohibitively expensive simulations, in terms of time and memory. As known from the theory of 
dynamical systems, seemingly complex high-dimensional dynamical systems can have low - dimensional 
global attractors, or low-dimensional invariant manifolds containing the physically relevant solutions of the 
system nimz]. In these cases, reduced order models (ROMs) exploiting the attractor or center manifold 
structure can dramatically reduce the time and memory needed to execute the corresponding full-order model 
(FOM). 

One of the most studied approaches to building reduced order models is based on Proper Orthogonal 
Decomposition (POD) [25l[T]. Normally, POD uses the (vector) solution of the FOM at selected time, space 
and parameter values, typically called ’’snapshots” ESlIl], to calculate a set of singular (orthonormal) vectors 
spanning the full space. In the context of POD, if the dimension of the FOM is n (a large number), for 
a given I < n, the span of the I vectors corresponding to the I largest singular values defines the reduced 
basis space (RBS). The ROM is defined and solved in the RBS and the prolongation of its solution back 
in n-dimensional space serves as an approximation to the solution of the FOM. The ultimate goal is to 
find an optimal (ideally small) I so that the approximate solution is sufficiently close to the solution of 
the FOM. Two typical problems are important for POD ROMs: (a) given a fixed error tolerance, find the 
the snapshots generating the smallest basis satisfying that tolerance; (b) given a fixed basis size, find the 
snapshots generating a basis that minimizes some error indicator. The specific choice of snapshots used 
affects greatly the size of the ROM basis and the ROM solution error. 

The availability of a priori error estimates or error bounds that include characteristics of the RBS (e.g., 
dimension, exact locations of the snapshots, type of snapshots, etc.) is important for understanding the 
origin of the approximation error and controlling it. A priori estimates could be used directly to determine 
the optimal reduced dimension of the model or to develop a posteriori error estimators to approximate it 
computationally. Although several papers deriving error bounds for POD ROMs have been published, e.g. 
0 nn [m na [m the bounds rarely contain information about the spacing between snapshots and none 
explicitly includes the exact locations of the snapshots. Including the latter information can be potentially 
useful for for informing snapshot selection and lower the computational resource requirements of POD. 
Notably, Kunisch and Volkwein m developed error bounds for POD ROMs for linear parabolic problems 
which are in the form of a sum of two terms: one depending on the singular values corresponding to 
the orthogonal to the RBS basis vectors and the other depending on the (uniform) distance between the 
time points where the snapshots were taken. This result was generalized and applied to other problems in 
subsequent publications, e.g. [IHE]. In the latter two publications the authors used a POD ROM based 
on two types of snapshots: a) collected from the solution and b) difference quotients (DQs) calculated from 
these snapshots. Although the DQs are in the span of the solution snapshots, their use was justified by a 
need to avoid the appearance of a blow-up term in the derivation of the error bound. 

A few other authors have used DQs or values of the time derivative as snapshots. In a variant of 
POD ROM, the Discrete Empirical Interpolation Method, applied to reduce nonlinear dynamical systems, 
Chaturantabut and Sorensen [Hi used (time-derivative) snapshots from the nonlinear part of the system. 
The incentive for using derivative snapshots in their method was to reduce the computational complexity of 
the POD ROM. However, so far, there is no clear answer whether using DQs or time derivative snapshots 
provides an advantage in terms of accuracy of approximation of the FOM solution. In a recent paper, Iliescu 
and Wang m posed the question whether there is advantage (and if so, what) to include the DQs in the 
POD calculation. The question was prompted by results of two groups of authors [HEs] who developed 
convergence analysis for POD ROMS that did not use DQs {i.e. providing evidence that using DQs is not 
necessary for avoiding blow-up of the POD error), as well as by conflicting numerical results on the beneht 
of using DQs in POD ROM formulations m- Iliescu and Wang defined a POD optimality criterion and 
explored the optimality, in several norms, of two POD methods - the first utilizing only solution snapshots, 
and the second - utilizing both solution and DQ snapshots. They found that the method with DQ snapshots 
was optimal in all norms while the one using only solution snapshots was optimal only in some of the norms. 

In this paper, we compare two POD methods for reducing large systems of ODEs. The first POD ROM 
method is based on snapshots from the solutions only; the second method uses snapshots from both the 
solutions and their time derivatives. The time derivative snapshots in general do not belong to the span of 
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the solution snapshots, therefore the number of singular vectors corresponding to non-zero singular values is 
normally twice larger {21) in the second case compared to the first (1). The additional computational cost 
to calculate the time derivative snapshots is negligible and thus, this is a cheap way to augment an existing 
basis; yet using derivative snapshots brings in additional information about the dynamics of the solutions. 
Thus, a method utilizing derivative information could produce more accurate ROM approximations which 
would be valuable when I is relatively small. These arguments provide an incentive for comparing these two 
methods. 

Our initial incentive, however, came from the expression for the error (j3.7p . where we noticed that if the 
vector f(x, t) was a snapshot, it would be in the span of the reduced basis vectors and the second term in 
the right-hand side of (13.71) would be zero. Intuitively, this property could make the error smaller compared 
to the case when this term is not zero. We explore this conjecture by deriving error bounds for the two 
methods. We have derived bounds in a form consisting of two terms: one depending on the largest neglected 
singular value and one depending on the distance between snapshots but not on the singular value. Assuming 
that the largest neglected singular value is small, we show that the bound of the 2-norm of the error has a 
second-order dependence on the length of the time intervals between snapshots when only solution snapshots 
are used, and a fourth order dependence on same length when also time derivative snapshots are used. 

To the best of our knowledge, this paper introduces two main new results. The first one is the actual form 
of the error bound which involves the time moments of the snapshots (via the time intervals between the 
snapshots) and the largest neglected singular value (instead of all neglected singular values). Error bounds 
and error estimates involving the time moments of the snapshots contain information which is potentially 
significant for rational snapshot selection. 

The second, and more significant, innovation is that the bounds allow for comparison between the two 
POD ROM methods (with and without derivative snapshots) and specifically suggest that if the first ne¬ 
glected singular value is sufficiently small, the method with time derivative information could be more 
accurate. Thus, these bounds give insights that we test numerically. Interestingly, we find that the behav¬ 
ior of the errors in numerical experiments with cases of the discretized FitzHugh-Nagumo system can be 
explained by the form of the bounds, although these bounds are not tight. 

The paper is organized as follows. Section 2 includes preliminary information and derivations. Section 3 
includes derivations of the error bounds for the two methods. Section 4 includes three numerical experiments 
illustrating the the validity of the insights obtained from the error bounds. 


2 Notations and relevant background 

Throughout the paper we will indicate vector- and matrix- valued variables by bold upper- and lower-case 
letters and scalar values by normal typesetting. We consider a dynamical system of the form 

X = f(x(t),<),x(0) = x°,t e [0,T], (2.1) 

where, respectively, boldface letters denote vectors, x(t) = {xi{t ),..., Xn{t))^ & M", f : —>■ M". 

Assumption 1. It is assumed that G C(R”) and ^ G C{0,T). Further smoothness assumption on f 
will be imposed in Proposition\^ 

Let y{t) be a solution of (12.ip . In the following theoretical treatment it is assumed that snapshots of the 
solution and its derivatives are collected at certain time points from y{t). We will explore the approximation 
error incurred by using a reduced order model (ROM) constructed from a POD basis using combinations of 
these snapshots. We first introduce some definitions and derive useful relationships. 


2.1 Using only y(ti) as snapshots. 

For the sake of self-consistency and setting up notations, we will revisit some well known definitions and 
results related to POD. 


We consider a set of time points, not necessarily equispaced, U G [0, T], i = 1,..., m < n where the latter 


inequality has been set mostly for clarity of presentation. We consider the matrix Y = 
with rank{Y) = mX. Then, 


y(ti):...:y(t™) 


G 
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(i) YY"^ is n X n matrix, Y^Y is an to x m matrix and rank(Y'Y'^) = ranfc(Y^Y) = mX < to. The 
squares of the to largest singular values of Y, {af )‘^,i = 1, to, are the eigenvalues of Y^Y, while the rest 
are equal to zero. 

(ii) Let uf = {uYi,U 2 i, G R" and vj = (u^, G M™ be orthonormal sets of eigen¬ 

vectors of YY^ and Y^Y. The vectors vt^ G R" and G are called left and right singular vectors 
of Y, respectively. Let and be the matrices whose columns consist of these eigenvectors, i.e. 
U^(U^)^ = I„, V^(V^)^ = Im- It is known that 


Y^uJ = alv. 




( 2 . 2 ) 


and 


■\r y Y Y ■ 1 

Yv, =a^ Ui ,1 = 1,...,TO, 


(2.3) 


where af > 0- 

Thus, YV^ = where is n x to matrix with diagonal entries aY > 0,i = and Os 

otherwise, and it is assumed that aY > aY > ... > ■ This leads to the well known singular decomposition 

of Y {e.g. [1] ): 

Y = U^S^(V^)^, (2.4) 

and the elementwise dyadic decomposition 


m 

yi*i) = '^<^I'>jYkuY. (2.5) 

fe=i 

(iii) In practice, once to snapshots are calculated, from them I < mX < m left singular vectors correspond¬ 
ing to the I greatest singular values of Y are used to construct the basis for the ROM. The rationale is based 
on the Schmidt-Eckart-Young-Mirsky theorem [T] which asserts that the truncated dyadic decomposition 


i 

i=l 


minimizes ||Y — X ||2 over all X such that rankX = 1. The value of I is typically chosen so that either the 
reduced basis approximates the snapshots with sufficient accuracy, or to satisfy constraints on the ROM 
(memory usage and computational time are both related to ^). From (12.51) the following 2-norm bound of 
the difference between the snapshot y(L) and its Z-truncated representation is derived: 








\ fc — /+1 fc — /+1 


A Y 

\ k^l+l 


\ Y i^k^Yy < 

\ k=l+l 


Y 

^l+l> 


\ k=l+l 


Y 


( 2 . 6 ) 


where (,) denotes the Euclidean dot product. 

In the sections that follow, we will use the following partitionings and notations. Given a matrix Y = 
U'^S^(V’^)^, we partition the columns of as follows. 






(2.7) 
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where = [ui:...:u;] is the reduced basis matrix , ] contains the truncated basis 

vectors with non-zero singular values and = [u^y_|_]^:...:u„] is the YY^ null space matrix: span(U^) = 
ker(YY^). Note that and may be empty matrices. 

Respectively, < mX are denoted as u^; if ^ < rrX, u;+i,...,u^y are denoted as = 

1, mX — 1; if rrX < n, are denoted as i = 1,n — mX. 

Since rank(YY^) = YY^U = 0„x(rt-mV) (the null n x {n — mX) matrix). Further, from (12.511 and 
the orthonormality of it follows Y^U^ = Omx(n-m)- 

2.2 Using both y(ti) and as snapshots. 

Consider a collection of 2m snapshots from both y(ti) and f (y(ti), U), i = 1, ..., m with 2m < n. (Again, the 
latter inequality is not a necessary condition, but rather imposed for clarity of exposition.) Let Z be the 
matrix made of these vectors. 


z = [y(ti):...: y(tm): f(y(ti), to): : f(y(tm), U)] G 


pn X 2m 


( 2 . 8 ) 


ZZ^ e and Z^Z e ]R 2 mx 2 m^ 

Let uf G R", i = 1,..., n and vf G i = 1,..., 2m be the orthonormal left and right column singular 


vectors of Z correspondingly, let = 
of singular values cti > ...ai-m- 
Then as in (I2.4I1 . 


uf :...: 


,V^ = 


rZ : 


: V 


2m 


and G matrix 


Z = U^S^(V^)^ , with rank(Z) = ' 


(2.9) 


and for I < , as in (1^ . 

i i 

\\yit^) < (jf+i and \\i{y{ti),U) - < af+i- (2.10) 

k=l k=l 

Similarly to the previous section, we denote by U^,U^,U^ the matrices corresponding to the eigen¬ 
vectors associated with the first I singular values, the I -I- l,...,m^ singular values and the nullspace of 

ITX . 

2.3 Projection estimates 

Let the column vectors G R", fc = 1,..., n be an orthonormal basis in R" consisting of the singular vectors 
of a snapshot matrix. In light of the above considerations, let ^ < m < n be given integers and let us define 
the matrices U = [u^] = [u^], for i = l,...,l; U = [u^] = [u^], for i = I + l,...,m; U = [ui] = [u^], for i = 
m + 1,..., n (one or both of the latter may be empty). 

Let us consider the n x n matrices (again, one or both of the latter may be empty) 

p = uu^, p = tru^, p = uu^. 

P is the projection onto spanjui} and P and P are projection matrices onto spanjuj} and spanjui}- 
We will be using of several occasions the following result which is easy to prove 

Proposition 2.1. P -|- P -|- P = I„ and ||P ||2 = ||P ||2 = ||P ||2 = 1- 


Let us consider the specific projection matrices derived from Y and Z, P^ = U^(U^)^, P^ = U^(U^)^, P^ 
U^(U^)^. The following estimates of the projections of the snapshots in Section [2TT] are derived as in (12.61) . 
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( 2 . 11 ) 


l|P’^y(i*)ll2 = = 

k^l 
mX —I 

||P’^y(i*)||2 = II (yl+ivJk+i^l\\2 = 


\ k=l 


k=l 


2'y-l 


. J2 (^^i+i^rk+i^ < ^i+i 

\ k^l 


and, 


l|P^y(iOll2 = o. 


( 2 . 12 ) 


(2.13) 


Similarly, if = TJ^('U^)^,P^ = tj^(U^)'^,P^ = (U^)^, the following estimates of the projec¬ 

tions of the snapshots from the solution and its time derivatives in Section [5] are derived as in (1^ . 


||P^y(i *)||2 < erf; 

||P^f(y(i.))||2<af; 

||P^y(i.)||2<af+i; 

||P^f(y(t.))||2<af+i; 
l|P^y(ir)ll2 = 0; 
l|P^f(y(i.))l |2 = 0, 

Having derived these estimates, in what follows, we will denote 

P + P = = I„ - P. 

Therefore, 

ll(P^)-"y(i 0 ll 2 = ll(P^ + P^)y(t.)l |2 = l|P^y(i,)ll 2 < for x = r,z. 


and 


P^)"-f(y(ir))ll 2 = ||(P^ -KP^)f(y(t .))||2 = ||P^f(y(t .))||2 < 


(2.14) 

(2.15) 

(2.16) 


3 Error bounds for the two ROM methods 

3.1 Derivation of an upper bound in the general case 

Let x(t) be any solution of (|2.1I) for some parameter set and initial conditions possibly different from those 
of y. Suppose that a number m of snapshots from the solution y have been calculated as described in the 
previous sections and a POD basis S K", i = 1, ■■■, n has been constructed. Since is a basis in R", x(t) 
can be decomposed as: 

x(t) = Ux(t)-f Ux(t)-f Ux(t), (3.1) 

where 

X = U^x S spanjui} S Rf x = U^x G spanjui} G R™”^ x = tj^x G spanjui} G R"”™ (3-2) 

are solutions of the dynamical system 

X = U’^f(i/x(t) -I- Ux(t) -I- Ux(t),t), x(0) = U^x(O) 

i = U^f(17x(t)-bUx(t)-GUx(<),t), x(0) = U^x(0) (3.3) 

X = U^f(17x(t) Ux(<) -t- Ux(<), t), x(0) = U^x(O). 

In the POD ROM approach the large n-dimensional ODE system (13.31) is replaced with an /-dimensional 
ODE system of the form 

z = U'^f(Uz(t), /), z(0) = U'^x(O), (3.4) 
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The prolongation of z in K”, Xp = Uz is the POD approximation of the solution x of the full system and 
satisfies the equation (as also stated in [H]): 


Xp(t) = Pf(Xp,t) 

Xp(0) = Px(0), 

while the exact solution of the full system satisfies 

i=(P + P^)f(x,t), 
x(0) = (P + P^)x(0),. 

Subtracting (13.5p from (13.61) . the equation for the error e{t) = x{t) — Xp{t) is 

e = P[f(x,t) - f(xp, t)] + P^f(x, t) 
e(0) = P^x(O). 


(3.5) 


(3.6) 


(3.7) 


We use the latter system to evaluate the error. As f is assumed to be continuously differentiable, 
f(x,t) — f(xp,t) = ^(x*(t), t)e(t), where for each t, x*(t) is a value defined by Taylor’s theorem such that 
is in the open interval with ends Xi{t) and xp ^{t). We denote 

A*.(t) = P^(x.(0,i) 

and note that it is dependent on the projection P also via Xp through x* (t). 

We note that Ap(t) is bounded on [0,T] and thus its 2-norm is bounded. Let Ap, be a Lipschitz constant 
such that 

||At(t)||2 < Ap,Vte [0,T]. (3.8) 

We next proceed to integrate system (13.71) (the integrals exist because of the boundedness): 

e(t) = P^x(O)-I- f P^i{x,s)ds+ f A*{s)e{s)ds. (3.9) 

Jo Jo 

Since f* P-'-f(x, s)ds = P-*“x(s)cis = P-*-(x(t) — x(0)), it is easily established that 

e(t) = P^x(t) -I- f PA*{s)e{s)ds. (3.10) 

Jo 

Applying 2-norms to both sides of (13.101) and the triangle inequality and noting that ||/^^ 0(s)ds||2 < 
la Il')’(®)ll2'^s ([ZI]) we obtain 


||e(f )||2 < ||PMi)ll 2 + r ||A*(s)e(s)|| 2 ds (3.11) 

Jo 

Taking into consideration Proposition 12.II and (13.8|) . the latter results finally in 

||e (<)||2 < ||P-^x(t)|| 2 -k Ap / ||e(s)|| 2 ds (3.12) 

Jo 

We will use this last inequality and a formulation of a Gronwall-type theorem from [2], Theorem 1.5.1, 
stating the following. 

Theorem 3.1. If g,A > 0 are real-valued continuous funetions on [0,T] and if the continuous function rj 
satisfies r]{t) < g{t) -I- A{s)r]{s)ds,t G [0,T], then 

r]{t) < g{t)-\- J A{s)g{s) exp(^J A{u) du^ ds,'dt G [0,T]. 
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Let 


g{t) = ||P*^x(t)|| 2 . 

Clearly g(t) is continuous. Applying Gronwall’s inequality, we get the estimate 

||e(t )||2 < 5(t) + J g{s)ds < ma,x^g{t). (3.13) 

The above is a bound for the difference between a solution of system m and a solution of the modified 
system ()3.5I1 . In deriving this bound, we have not used so far how the projection matrix P was constructed. 
Obviously, the error will be zero if the solution x(t) is orthogonal to P'^, i.e. if x(t) is in the span of the 
basis vectors comprising U. Clearly, if the snapshots used to define P are chosen so that x(t) is P^x(t) 0 
over the whole interval [0,T] {e.g. the solution snapshots are sufficiently dense and the solution does not 
change rapidly in time) it is obvious that the error will be small. It also makes sense to consider snapshots 
from the time derivatives, or, possibly, other characteristics of the solution that contribute to reducing the 
term g{t). We explore this idea in the next sections. 


3.2 Error bounds for the approximation of the solntion from which the snap¬ 
shots were collected. 

One could use the reduced basis derived from a particular solution of the FOM (12.11) . It is important to 
derive error bounds in these cases; however, here we start with deriving an a priori bound of the error made 
when approximating the solution y from which the snapshots were collected. We will derive error bounds 
for the two ROMs: a) when only snapshots from y are used; b) when snapshots from both f and y are used. 
To this end, we need to make an estimate for g{t) in (I3.13p . 

Let tj, i = 1,..., m < n,ti = 0,tm = T he the time points at which the 2m snapshots y(ti) and f (y(ti), U) 
were calculated. Denote = ti+i — 


3.2.1 Method 1 - Using only solution snapshots 


Let the snapshot matrix be Y = [y(ti) :...: which, given I < m, generates the projection matrix P^ 

and let the solution of (13.51) be ypvit). 

Let us take some value of t G [0, T] and let it belong to the interval ti+i] for some i. Applying Lagrange 
interpolation, we get 


y{t) = y{u)^ —+ y(^i+i)7r-—where 


{ti — ti+i) 


{ti+i — ti) 


R(yW) =-;;-- 


(3.14) 


dt 


where C G (tidi+i) and where it is assumed that f has first derivatives in all variables [53]. Because of 
Assumption jT] is bounded on [ti,ti+i], i.e. there exists a constant 4'^, such that 


, df{y{t),t) 

dt 


I 2 < 4**,^ G [U,ti+i]. 


Multiplying the above expression for y(t) by (P^)-^, we get 


(P'^)-"y(t) = (P'^)^y(t,) 


t — t. 


i+l 


+ (P ) yiU+i)' 


t - ti 


ti ti.^1 ti 

Ai 


- ti)it - U+i) di{t) 
w ) -^- —\t=C 


dt 


Since maxtgfj.^t.^j] \{t-ti){t- ti+i)\ = ^, it follows 

||(P^)^y(i )||2 < ||(P^)-^y(t0ll2 + ||(P^)-^y(t,+i )||2 ■ 

Using ()2.15|) we get from 1)3.17^ : 


^1 di{t) I 


(p"r 


dt 


|t=Cl|2 


(3.15) 


(3.16) 


(3.17) 


g{t) = ||(P^)V( 0 ll 2 < + \\{P^ )"||2^4/,. 


Y 




(3.18) 



















Using (13.131) we get: 


WII 2 = l|yW -ypr(t)||2 < 


2a, 


f+l 


A2 


^pY t 


(3.19) 


where aY^i is the {I + l)-th singular value of the snapshot matrix Y and and Apr are defined via p.l5l) 
and dSH) (with P = P^). 


Thus, we obtain the following 

Proposition 3.2. Let f satisfy Assumption 1 and let y(t) be a solution of i2.1\) . Let tj,j = 1, ...,m < n 

be a set of points in [0,T] and let Y be a matrix of solution snapshots, Y = [y(ti):... :y(tm)]- Let 
be the matrix of the truncated set of the first I < m singular vectors u^,fc = l,...,l of Y and P^ be the 
corresponding projection matrix; P^ = U^(U^)^. Let yp^ft) be a vector function whose projection in 
span{u^,fc = 1, ...,1} solves the ROM Then the 2-norm of the error (t) = y{t) — yprit) satisfies 

the bound iS.ltA) . where — ti and t € [U,<i+i]. 


3.2.2 Method 2 - Using snapshots of the solution and the time derivatives 

We use similar logic to the one in the previous section. 

Let the snapshot matrix be ^ = [y(to) ^ ^ y(tm) ■ f(y(io), ^o), ^ ■^{y{tm),tm)], which, given I, generates 

the projection matrix P^ and let the solution of (13.51) be ypz(t). 

Let t G [0,T] belong to the interval [U,U+i] for some i. We now apply Hermite interpolation to define a 
vector polynomial p*(t) such that p*(tfc) = y(tfe), ^\t=tk = f(y(4),tfc),fc = + 1. 


P*(^) =yiti) 


2{t-U) 1 {t-U+if 

(^2+1 ^i) ^ y^i ^2+1) 


1 + 


2(i-t,+i)] {t-Uf 


f(y(U),U) *'^ ^, 2 ^^ + f(y(U+i),U+i) 

[ti - ti+iY 


[ti -ti+i)i [ti -ti+i)2 
[t - tiY[t - tj+i) 

[ti - ti+l)2 


(3.20) 


and 


y(t) = p*(t) + ^[t - tif[t - 


(3.21) 


where 9 G (ti,U+i) and where it is assumed that f is three times differentiable in y and in t on [0,T]. 
Assumption 2. We further assume that || II 2 is continuous and thus, bounded on [0,T]. 

Thus, for each interval \ti,ti^i\, there exists a constant such that 

,,dfi[t) 


dt^ 


-II2 [o,r] 


We apply the following bounds 


(3.22) 


max |2(t - U)[t - U+if\ = max |2(t - tif[t - ti+i)\ = —Af; 
te[ti,ti+i] tG[ti,ti+i] 27 

A2 

max [t — tiY= max [t — ti+if = —Y, 

A^ 

max \[t- tif[t - ti+if\ = 


and (j2.15112.161) in (|3.20l) and (|3.21|) to get the bound 


Af 


g[t) =||(P^)^y(t )||2 < [2a4i(- + - + -A,) + 




(3.23) 


(3.24) 


9 































Using this upper bound, we obtain from (13.131) : 


l|e^(t)||2 = l|yW -ypzWIb < 


-'/+ 1 ' 


,59 

M 


4 

_A4 + 

27 ’ 384 


(3.25) 


where is the I + 1-th g of the snapshot matrix Z and <i>i and Apz are defined via (13.221) and (13.81) . Apz 
is the bound for the Jacobian as defined in section lO and depending on the specific projection . 

Thus, we obtain the following 

Proposition 3.3. Let f satisfy Assumptions 1, 2 and let y(t) be a solution of \2.1\) . Let tj,j = 0, ...,to < n 

he a set of points in [0, T] and let Z = [y(to) ■ • y(lm) '■ f(y(^o)) ^o)) • ••• • f (y(im)j ^m)], *-e., Z is a matrix of 

solution and time derivative snapshots. Let \J^ be the matrix of the truncated set of the first I < 2m singular 
vectors , k = 1, ..., I of Z and be the corresponding projection matrix; P^ = U'^(U^)^. Let ypz(t) be 
a vector function whose projection in span{u^,fc = solves the ROM Then the 2-norm of the 

error (t) = y{t) — ypzit) satisfies the bound iS. 25\) . where — U and t € [ti,ti+i]. 


3.3 Remarks on the error bounds 

(a) To the best of our knowledge, the bounds (13.191) and (13.251) are the first derived bounds that include the 
time points at which the snapshots were taken. So far published bounds 0 0 HU [13 m ng, if including 
time information at all, assume equidistant snapshots and include the time step. Our bounds do not assume 
any specific distribution of the snapshot times. For equidistant snapshots A^ = A, these bounds also contain 
local information by including derivative bounds in the intervals Including local information may 

be helpful for rational snapshot selection via error estimates. 


(b) We note that the derived bounds are still far from being exact or being fully informative. For a linear 
system x = Ax, 


e{t) = P^y(O) -k / e-^"^P^Ae"^'*P^y(0)ds 


(3.26) 


i.e. the error depends on the eigenvalues of PA, which is neither reflected in the bounds derived, nor, in fact 
in any bounds or estimates published. Yet, the derived bounds are valuable because they provide insight 
about the relationship between the basis truncation and the location of snapshots. 


(c) For linear systems of the form x = Ax + b, an explicit bound, not depending on P, can be found. Indeed: 


||A^(t)||2 = ||PA||2 = ai(PA) < IIAlb = ai(A), (3.27) 

where cri(A) is the largest singular value of A, i.e. we can assume that Ap = cri(A). Also, let 9i = 
I|y( 0 ll 2 - Then 

Ti = max ||Ay(t )||2 < ||A|| 26 (i = ai{A)0i, 

= max ||A^y (<)||2 < ||A||f6»i = (cri(A))^6»i 


and from (13.191) and (13.25^ : 


l|e^Wll2< 


2 a, 


i+i 


A 2 


„cri(A)t 


(3.29) 


and 


l|e^Wll2< 


^i+i 


118 4 

108 ^ 


A,) 


A4 

(o-i(A))3^6», 


„CTi(A)t 


(3.30) 


We next try to compare the bounds for the two ROMs considered. Above, the only quantities that 
depend on the method used are and crffi. Obviously, if the dimension of the RBS is taken to be equal 
to the number of snapshots used (i.e. I = m in the first case and I = 2m in the second case and then 
= 0), the reduced model by Method 2 will have twice larger dimension than from Method 1, 
but, for small A^, its error might be much smaller (possibly with two orders in A^ as predicted by the error 
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bounds). Thus, for small ai+i and A^, Method 2 may lead to a significantly smaller error. Note that the 
terms ’’possibly” and ’’may” are used because these are upper bounds and not exact estimates. We check 
these predictions in the next section. 

In the case where Ci+i is not sufficiently small, using both solution and derivative snapshots yields an 
error bound that is hrst order in A^, whereas using solution snapshots only yields an error bound that is 
zero-th order in A^. So improvements in the error for the solution snapshot - only case come primarily from 
increasing the dimension of the RBS, whereas the error in the the time derivative and solution snapshot case 
decreases when snapshots are sampled more frequently. 

These derivations demonstrate the trade-off between error and dimension of the RBS for the two methods 
and call for investigating the comparison between the distributions of the singular values for the two methods. 
For example, for a fixed I and equidistant snapshots, when A^ = A is decreased, (7^^^ and crf^-^ may increase 
(since adding more snapshots means adding more singular values). Therefore, there may exist an optimal 
value of A such that further decrease of the distance between snapshots would not decrease the error as it 
will be dominated by the error caused by the truncation of the snapshot-generated basis. We demonstrate 
this phenomenon in the next section. 

Now let us look at the case when the dimension I of the RBS is fixed for both methods a.t I = m. Then 
o'Y+i = 0 but > 0. However, <jY\.i may be insignificant compared to the rest of the error. If the bounds 
were good estimates of the error, the error of Method 2 would be smaller than the error of Method 1 in such 
cases. This possibility is explored and demonstrated on examples presented in the next section. 

Arguments similar to the ones presented above for linear systems hold also for nonlinear systems for 
which the existence of a constant A not dependent on P, such that A> Ap can be proved. 

(d) As pointed above, investigating the distribution of the singular values and how it changes with 
decreasing A is important. It is not clear how dependent on the particular problem this distribution is. For 
the purpose of comparing methods using snapshots containing additional information about the problem 
(such as derivatives), it is important to understand how adding this information may change the distribution 
of singular values. Some examples demonstrating these differences are considered in the next section. 

4 Numerical experiments 

To validate the above bounds and to compare numerically the error from the two methods, we performed 
numerical experiments with systems of ODEs derived from a method-of-lines discretization of the FitzHugh- 
Nagumo system with diffusion (FHND) [13]. This system has been used as test problem in various studies 
of model reduction methods, e.g. |5]. To provide some background, the FHND system is an approximation 
of the Hodgkin-Huxley system of equations, designed to describe the propagation and dynamics of an action 
potential (difference in external and intracellular voltage) generated along the nerve axon. Since the dynamic 
behavior of the solutions to the FHN system is very sensitive to changes in some of the parameter values 
[14] . and the solutions are characterized by a combination of fast and slow dynamics, it is often used as a 
test case for the accuracy of numerical methods. In the FHND system, V(x,t) is the membrane potential, 
ie., the difference between the extracellular and intracellular potentials and w(t) is a ’’recovery variable”. 

Different texts (e.^., HZ], HO] , [SI) consider different forms of FHND. In general, the ID version of the 
system is a system of two reaction-diffusion equations 


vt = DiVxx + fi{v,w) 
Wt = D2 Wxx + f2{v,w) 


(4.1) 


where x € [0, A], representing a one dimensional axon with length X. In most texts, fi is cubic in v and 
linear in w: /i = A u(l — v){v — a) — w and f 2 = cv — bw. 

The initial conditions are set up so that initially the nerve membrane is at equilibrium: 


w(a:, 0) = 0; w(a::, 0) = 0. (4.2) 

Different boundary conditions (BC), depending on the problem, can be considered as outlined in [13] . 
The Neumann BCs Vx{0,t) = —Io{t),Vx{X,t) = 0 correspond to applying current Iq at the ’’left” (at 0) end 
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of the axon (where r is a constant, depending on the internal and external resistances) and ’’sealing” the 
axon at the other end (no current). 

In this paper, we solve numerically equations (14.11) together with the following boundary conditions for 
w{x, t)\ 

w{0,t) = woit) 
w{X,T) = wx{t), 

and the BC for v{x,t): 


Vx{0,t) = -loit) 
VxiX,t) = 


(4.4) 


Above, WQ(t),wx(t), lQ{t), Ix{t) are input functions. 

Equations (10-01-01) were semidiscretized using the method of lines with finite differences. We 
denote: n = 2(L + 1); Ax = X/L;xj = j Ax, j = 0,L-,Vj(t) = v{xj,t),j = 0, ...,L 
Using the approximations 


.( 0 ,t) 


V2it) - Vl{t) 

Ax 


- VxiO,t) 


Vx{X,t) - 


Vxx 


Ax 

VL - VL-1 

Ax 


7 ^XX\^J 


{Xj,t) 


Vj-i + Vj+i — 2vj . 


(Ax )2 


,j > 0; 


(4.5) 


Ax 


we derive discretizations by the method of lines for the equations (14.11) - (14.3114.41) as follows. 


^ = (^)2 ['*^2 -vi + Ax/o(t)] + fi{vo,wo), 

dvj Di 

~dt ^ + /i(uj,u;j),j = 1,...,L- 1 

^ = (^^)2 ['^^-2 “ - Ax/i(t)] + fl{vL,WL) 

dwj D 2 

~dr ^ (Ax)2 ^ j = 1,..., L - 1 

where the initial conditions are Vi{0) = 0, Wi{0) = 0 and wo,wl are defined in (14.31) . 


(4.6) 


4.1 Experiments 

The two POD ROM methods described above were implemented in a Matlab code. Specifically, we solved 
equations (14.31) . (14.4L (14.61) and selected m equally spaced on the time interval [0,r] snapshots from the 
solution. The ODE solutions were obtained numerically using Matlab’s routine odes 15s, where the absolute 
and relative tolerances were sufficiently small (usually equal to 10“^^, 10“^"'’) to ensure stable performance 
of the integrator. The right hand side of (14.61) was calculated using the selected snapshot vectors to obtain 
the time derivative snapshots. The snapshot matrices Y and Z were formed and their SVD calculated using 
Matlab’s svd routine. The dimension I of the ROM was either predefined or chosen so that ai+i < e < ai 
for a predefined e. After determining the value of I, the respective projection matrices for the two methods, 
and were calculated as described. The solutions of (13.5|) with P = P^ and P = P^ were calculated 
respectively for the two models. The distribution of the singular values in the two methods for different 
snapshot selections was also calculated. 

To be more specific, let us denote by = (v^,w^) the solution obtained by using only solution 
snapshots and by = (v^, w^) the solution obtained by using both solution and derivative snapshots. We 
calculate ||e^(t )||2 = ||y(t) -y'^(t )||2 and ||e^(t )||2 = ||y(t) -y^(t)|| 2 - 
All experiments were done with X = 10, wq = wi = 0. 
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4.1.1 Experiment A 


The two ROM methods were explored for a linear homogeneous system, obtained from (14.61) with A = 0, of 
402 (ie., L = 200) equations with constant coefficients Aa; = 0.05, £>1 = 15, D 2 = 10,/r = 10,7 = 5 and 
Io{t) = = 5. ODE solutions were computed in Matlab using odel5s with both absolute and relative 

tolerances set to 10“^^. The dynamics of the system can be predicted theoretically. It has no equilibria, the 
eigenvalues are nonpositive and there are 2 zero eigenvalues. Due to the zero initial conditions, the solutions 
Vj,Wj behave asymptotically like where 5 > 0 is a constant, so the solutions eventually converge to 0 . 

The system was solved on the time interval [0,0.5]. The calculated FOM solution is shown on Figure SI 
(Supplement). 

For both methods, three sets of equidistant snapshots were collected, with A = 0.01 (50 snapshots), 
0.005 (100 snapshots), 0.0025 (200 snapshots) and the errors ||e^(t)|| 2 , ||e^(f )||2 were calculated. The plots 
in Figure 1 demonstrate the size of these errors when A is reduced while keeping the first neglected singular 
value CTi+i constant, to compare with the predicted behavior of ||e^(t )||2 and ||e^(t)|| 2 . 

The top two plots in Figure 1 correspond to ai+i < 10“^®. If the bounds (13.191) and (13.251) were 
exact predictors of the error, the distances between the curves for these cases would be roughly equal to 

logic lleXWIb - logic l|eX/ 2 ( 0 ll 2 = logic ||ef^‘;|(|)|| ~ logic = logic( 2 ^) - 0-602 and log^o ||ei(t )||2 - 


logio I|oa/ 2 (^)II 2 ~ logic(2^) ~ 1.204, respectively (in the latter the sub-index A was added to denote the 
value of the spacing between snapshots used in the calculations). The actual distances appear to be twice 
as big, indicating that the error is probably of higher order in A than the error bounds. 

For larger ai+i the distances between the curves decrease below these values (0.602 and 1.204). This 
behavior would be expected if the bounds were exact estimates of the error. Indeed, if this was the case 


logic INa Wlb-logic l|eX/ 2 WI |2 = logic l]kf 


A/2 


(i)ll 


logj^Q /4 is decreasing function of <7i+i. 


This behavior is demonstrated on the middle two (crz+i < e = 10“®) and bottom two (crz+i < e = 10“^) 
plots. Increasing the value of ai+i in general decreases the distance between the curves corresponding to 
different decreasing values A and increases the error in all three cases since it is dominated by the value of 
CT;+i. Note that the error corresponding to smaller A is not always smaller (middle and bottom right plots) 
when the value of di+i is significant. 

Often in practice, the dimension of the ROM is predefined. Therefore, it is interesting to compare the 
performance of the two methods with fixed ROM basis dimension. We present the results of an experiment 
with the same linear system where we compare the error from the two methods with the three different time 
steps and with a fixed dimension of the RBS in Figure 2. For the 402-variable FOM, equidistant snapshots in 
time of both the FOM solution and its time derivative were collected at spacings of A = 0.01,0.0025,0.005. 
After calculating the respective singular vectors, the first I = 5,10,..., 50 vectors were used to calculate the 
ROM solution. Presented are 6 of the calculations, with ROM basis dimensions 5, 10, 15, 20, 35, 50 to 
illustrate the observed tendency (Figure S2). The plots in red in Figure 2 correspond to the three ROMs 
via Method 1 and the plots in blue correspond correspond to ROMs calculated by Method 2. 

To get understanding of the plots presented on Figure 2, we investigate the distributions of the base-10 
logarithm of the singular values in the two ROMs (Figure S3). 

For I = 5, all three types of snapshot spacings in Method 2 produce extremely large error (> 10®), while 
the error produced by Method 1 is of size (10“^) comparable with the amplitude of the solution (0 — 0.5). 
Concomitantly, the first neglected singular value, aj, for all 3 snapshot selections in Method 1 is less than 
than 1, while , for all 3 snapshot selections in Method 2 is at least 10 times larger (Figure S2). In this 
case, for both methods, decreasing the distance between snapshots has insignificant effect on the value of 
the error, which is evidently dominated by the neglected singular value. 

When the ROM basis dimension is increased to I = 10, ||e^(t )||2 and ||e^(t )||2 are comparable. For 
I > 15, the error produced by Method 2 is considerably smaller than from Method 1. Having in mind the 
derived bounds, we explain this behavior with the steep decrease of and and consequent dominance 
of the 0(A®) and 0{A'^) in Methods 1 and 2, respectively. At Z = 15 the neglected singular value is less than 
10“^, at Z = 25 it is less than 10“®, and at Z = 35 it is less than 10“^® for all between-snapshot distances 
for both Method 1 and Method 2 (Figure S2). For this problem, for ROM dimension higher than 15, the 
error is mostly due to the terms that do not depend on the truncation of the basis, but depend on the 
error contributed by the spacing of the snapshots; therefore the method with higher order (0(A‘^)) produces 
smaller error. 
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Decimal log of errors, N1 = 200, N2= 100. N3= 50,n=402 Decimal log of errors, N1= 200, N2= 100, N3= 50,n=402 

BCy0t=-r0;y1t=-r1, r0=1.00 r1= 5.00, x= 10.000 BC y0f=-t0; y1t=-r1, r0=1.00 r1= 5.00, x= 10.000 



Figure 1: Experiment A. Error from the two methods at three different values of A = 0.01,0.0025,0.005 and 
different cutoff (e) values. The x-axis is time t and the y-axis is logiQ(||e^(t)|| 2 ) (left) and logiQ(||e'^(t)|| 2 ) 
(right). Circles correspond to At = 0.01, dots - to At = 0.005, and crosses to At = 0.0025. The plots on the 
right correspond to error from Method 2 (solution and derivative snapshots) and plots on the left correspond 
to error from Method 1 (no derivative snapshots). 
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Decimal log of errors, N1 = 200, N2= 100, N3= 50, n= 402 Decimal log of errors, N1 = 200, N2= 100, N3= 50, n= 402 

BC y0f=-r0; yn=-r1, r0=1,00 r1= 5.00, x= 10.000 BC y0f=-t0; yn=-r1, r0=1.00 r1= 5.00, x= 10.000 

a=0.00 D1=15.00 D2= 10.00 lambda= 0.00 gamma= 5.00 mu= 10.00 10= 5, 11= 5, 12= 5 a= 0.00 D1=15.00 02=10.00 lambcla=0.00 gamma= 5.00 mu=10.00 10= 10, 11= 10, 12= 10 



Figure 2: Experiment A. Error from the two methods at three different values of A = 0.01,0.0025,0.005 and 
different fixed RBS dimensions {I = 5,10, 35, 50). The x-axis is time t and the y-axis is logiQ(||e^(t)|| 2 ) (left) 
and logjQ(||e^(t)|| 2 ) (right). Circles correspond to A = 0.01, dots - to A = 0.005 and crosses to A = 0.0025. 
Red - Method 2, blue - Method 1. The y-axis is the decimal logarithm of the error. 
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We note again that we have only derived upper bounds and not exact estimates of the error. The phe¬ 
nomenon demonstrated above may not necessarily be valid in all cases. However, Experiment A demonstrates 
the potential for achieving better accuracy of approximation when using Method 2 instead of Method 1 for 
relatively low-dimensional ROMs. 

4.1.2 Experiment B 

In this experiment, the two ROM methods were explored for the nonlinear system (j4.6ll with A = 2, of 402 
{i.e. L = 200) equations with constant coefficients Ax = 0.05, Di = 5,D2 = = 1,7 = 5, a = 0.1 and 

Io{t) = 1.5(sin<)^, Ji(t) = 0.5(sint)^. Again, ODE solutions were computed in Matlab using odel5s with 
both absolute and relative tolerances set to 10“^'^. 

For the convenience of the reader, some of the text is identical with the previous section. 

The system was solved on the time interval [0,2]. The calculated FOM solution is shown on Figure S4 
(Supplement). The figure shows the plots of the 402 variables yj{t) where time t is on the x-axis. 

As in the previous example, for both methods, three sets of equidistant snapshots were collected, with 
A = 0.04,0.02,0.01 and the errors from Method 1 and Method 2, ||e^(t)|| 2 , ||e^(t)|| 2 , were calculated. The 
plots on Figure 3 demonstrate the size of these errors when A is reduced, while keeping the first neglected 
singular value ai+i below a constant threshold, to compare with the predicted, by the derived bounds, 
behavior of ||e^(t )||2 and ||e^(<)|| 2 . 

The top two plots in Figure 3 correspond to cr;+i < 10“^^. Even though the problem is nonlinear, the 
distances between the error curves suggest, similarly to the previous numerical experiment, that the error is 
of higher order in A than the error bounds indicate. The error incurred by applying a ROM with derivative 
snapshots is again at least one order smaller than the error incurred by Method 1. 

The middle and bottom two plots demonstrate the dependence of the error on the size of the largest 
neglected singular value (cTi+i < e = 10“^ and (T;+i < £ = 10“^) respectively. Similarly to the previous 
numerical experiment (A), increasing the value of (T;+i generally decreases the distance between the curves 
corresponding to different decreasing values of A and increases the error in all three cases, since it is dom¬ 
inated by the value of (T;+i which corresponds to the behavior of the error predicted by the error bounds. 
This behavior is clearly seen on the bottom plots of Figure 3, where the plots indicate that reducing A has 
miniscule (right plot. Method 2) or zero (left plot, method 1) effect on the magnitude of the error. Note 
again that the error corresponding to smaller A is not always smaller (middle and bottom right plots) when 
the value of cr/+i is significant. 

Next we present the results of an experiment with the same nonlinear system where we compare the 
error from the two methods with the three different snapshot spacings where the dimension of the RBS is 
fixed (Figure 4). For the 402-variabie FOM, equidistant snapshots in time of both the FOM solution and 
its time derivative were collected at spacings of A = 0.01,0.02 and 0.04. Dimensions of 5,10,15,..., 50 were 
considered. Presented are 4 of the calculations, with ROM space dimensions =5, 20, 25, 50, to illustrate the 
observed tendency. For a very low ROM dimension, e.g. I = 5, Method 1 and Method 2 produce similar 
errors (ss 10“^) and for each method decreasing A does not have any effect on the error (Figure 4). At 
Z = 5, (Te ~ 10“^ for both methods and all three values of A (Figure S4, Supplement) and the magnitude 
of the error is determined by the basis truncation as suggested by our error bounds. We further focus on 
the case with A = 0.04 (50 collected snapshots, plots with circles). When the ROM basis has dimension 
I = 20 ||e^(t )||2 becomes larger than ||e'^(t )||2 and remains so for all fixed RBS dimensions 25, . . ., 50. 
The same tendency of the error from Method 2 to become smaller than that of Method 1 when increasing 
the RBS dimension is observed in the case with A = 0.02 (100 snapshots). More precisely, for the case when 
100 solution snapshots were collected, at I = 25, ||e^(t )||2 becomes larger than ||e'^(t )||2 and remains so for 
the larger RBS dimensions 30, 35, . . .,50 (dotted plots). This tendency is likely true for the third case 
A = 0.01 (200 snapshots) as well (calculations with RBS dimensions over 50 were not done). 

The latter again demonstrates the potential for achieving better accuracy of approximation when using 
Method 2 compared to Method 1 for relatively low - dimensional ROMs. 

Similarly to the previous example, plotting the distributions of the singular values helps to get insight 
about the behavior of the error, observed on Figure 4. The plot of the distributions of the singular values 
in the two ROMs is presented on Figure S6 (Supplement). It is observed that for a given A, there is a 
dimension I such that cr/+i is comparable to the magnitude of the error produced by both methods. Then 
Method 2 produces less error than Method 1 for ROM dimensions greater than 1. 
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Decimal log of errors. Ml = 200, N2= 100, N3= 50, n= 402 Decimal log of errors, N1= 200, N2= 100, N3= 50, n= 402 

PeriodicBCy0t=-r0'(sin{t)f, y1t=-rr{sln(t))^, r0=1.50 r1=0,50, x= 10.000 Periodic BC y0t=-r0'(sln(f)f, y1t=-rr(sin(t))^ r0=1,50 r1=0.50, x= 10,000 

eps= 1e-15 a=0.10 D1=5.00 D2=1.00 lambda=2.00 gamma=5.00 mu=1.00 10= 95, 10= 57, 10= 41 eps= 1e-15 a=0,10 D1=5.00 D2=1.00 lambda=2.00 gamma= 5.00 mu= 1.00 11=224, 11=139, 11= 85 



Figure 3: Experiment B. Error from the two methods at three different values of At = 0.01,0.02,0.04 and 
different cutoff (e) values. The x-axis is time t and the y-axis is logiQ(||e^(t)|| 2 ) (left) and logiQ(||e'^(t)|| 2 ) 
(right). Circles correspond to A = 0.04, dots - to A = 0.02 and crosses to A = 0.01. The plots on the right 
correspond to error from Method 2 (solution and derivative snapshots) and plots on the left correspond to 
error from Method 1 (no derivative snapshots). 
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Decimal log of errors, N1= 200, N2= 100, N3= 50, n= 402 Decimal log of errors, N1= 200, N2= 100, N3= 50, n= 402 

Periodic BC y0t=-r0*{sin(t))^, y1t=-r1*{sin(t))®, r0=1.50 r1= 0.50, x= 10.000 Periodic BC y0t=-r0*(sin(t)f, y1t=-r1*{sin(t))®, r0=1.50 r1= 0.50, x= 10.000 

a=0.10 D1=5.00 D2=1.00 lambda= 2.00 gamma= 5.00 mu=1.00 i0= 5, 11= 5, 12= 5 a=0.10 D1=5.00 D2=1.00 lambda=2.00 gamma=5.00 mu=1.00 10= 20, 11= 20, 12= 20 



Figure 4: Experiment B. Error from the two methods at three different values of At = 0.04,0.02,0.01 and 
different fixed RBS dimensions {I = 5, 20, 25, 50 and for the problem in Experiment B. Circles correspond to 
A = 0.04, dots - to A = 0.02 and crosses to A = 0.01. Blue - Method 1, red - Method 2. The y-axis is the 
decimal logarithm of the error. 

For A = 0.04 and I > 25, ct;+i < 10“® for for both Method 1 and Method 2 ROMs and the maximum of 
error is of similar order (10“®'^ — lO”"^'^. As we saw, at I > 20, Method 2 yields more accurate approximation 
of the FOM than Method 1. As the bounds we derived suggest, this phenomenon is a trade-off between the 
value of cr/+i which at some point (here I = 20) becomes sufficiently small so that the error is dominated by 
the terms containing A^ and A^ in both methods. The latter trade-off leads to Method 2 becoming more 
accurate than Method 1, because of its higher order of approximation O(A^). 

4.1.3 Experiment C 

In this experiment we explore what happens with the error if the distance between the times at which the 
snapshots were collected is not small. We consider the same nonlinear system of equations as in experiment B, 
which was, however, solved numerically on the time interval [0,20]. For the ROMs 40, 20 and 10 equidistant 
snapshots were collected with A = 0.5,1, 2 respectively. The solution of the system is shown on Figure S5 
(Supplement). 

We consider the same type of scenarios as in the previous examples. Having in mind the bounds we 
derived, we still expect to see decrease in the error when cr;+i is very small. This is indeed demonstrated 
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on Figure 5, top plots (a-;+i < 10“^®). Again, as in the previous examples, the difference between the plots 
is larger than expected. Similarly to the previous experiments as well, for larger the errors become of 
similar magnitude when A is decreased, indicating that the value of (J;+i dominates the error. Interestingly, 
when t grows, the error decreases initially and then stabilizes around some value. For this specific problem, 
the error produced by Method 1 is larger than the error from Method 2, even for A > 1. The error is 
remarkably small given the relatively small number of snapshots taken. 

Figure 6 presents the results of an experiment with the same nonlinear system where we compare the 
error from the two methods with the three different time steps and where the dimension of the RBS is fixed. 
Presented are plots of the error for 8 cases where the RBS dimensions are held fixed at 5, 10, 15, 20, 25, 30, 
35 and 40. When the dimension of the RBS is 5 or 10, the magnitude of the error is similar for both methods 
and all distances between the snapshots (top two plots on Figure 6) and also comparable to the magnitude 
of the first neglected singular value respectively). For larger RBS dimensions, the error produced 

by Method 2 becomes smaller than that produced by Method 1. Specifically, since one cannot construct a 
POD ROM of higher dimension than the number of snapshots collected, the approximation from Method 2 
becomes better than that from Method 1 when increasing the the RBS dimension above 10 (circle symbol 
plots). 


Decimal log of errors, N1= 40, N2= 20, N3= 10, n= 402 
PeriodicBCy0t=-rt)*(sin(t)f,y1t=-r1*(sin{t))^ r0=1.50 r1=0.50, x= 10.000 
eps= 1e-15 a=0.10 D1=5.00 D2=1.00 lambda=2.00 gamma= 5.00 mu= 1.00 10= 40, 10= 20, 10= 10 



Decimal log of errors, N1= 40, N2= 20, N3= 10, n=402 
Periodic BC y0f=-t0*(sln(t)f, y1f=-rr(sln(t))^, r0= 1.50 r1=0.50, x= 10.000 
eps= 1e-15 a=0.10 D1=5.00 D2=1.00 lambda=2.00 gamma= 5.00 mu= 1.00 11= 80, 11= 40, 11= 20 



Decimal log of errors, N1= 40, N2= 20, N3= 10, n= 402 
Periodic BCy0t=-r0*{sin(t)f, y1t=-rr{sin{t))^ r0=1.50 r1=0.50, x= 10.000 
a=0.10 D1=5.00 D2=1.00 lambda=2.00 gamma=5.00 mu=1.00 10= 17, I0= 13, 10= 10 


Decimal log of errors, N1= 40, N2= 20, N3= 10, n= 402 
Periodic BCy0t=-r0*{sin(f)f, y1t=-rr{sin{t))^, r0=1.50 r1=0.50, x= 10.000 
eps= 0.0001 a=0.10 D1=5.00 D2=1.00 lambda=2.00 gamma=5.00 mu=1.00 11= 20, li 



Figure 5: Experiment C. Error from the two methods at three different values of At = 0.5,1, 2 and different 
cutoff (e) values. Circles correspond to A = 0.5, dots - to A = 1 and crosses to A = 2. Left plot - Method 
1; right - Method2. 
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Decimal log of errors, N1= 40, N2= 20, N3= 10,n=402 Decimal log of errors, N1= 40, N2= 20, N3= 10,n=402 

Periodic BC yOt=-rO*(sln(t)f, y1t=-r1*(sin{t))® r0=1.50 r1=0.50, x=10.000 Periodic BC yOt=-rO*(sin(t)f, y1t=-rr(sin{t))® r0=1.50 r1=0.50, x=10.000 

a=0.10 D1=5,00 D2=1.00 iambda= 2,00 gamma=5.00 mu=1,00 i0= 5, i1= 5, i2= 5 a=0.10 D1=5.00 D2=1,00 iambda=2.00 gamma= 5,00 mu=1.00 i0= 10, i1= 10, i2= 10 




Figure 6: Experiment C. Error from the two methods at three different values of A = 0.5,1, 2 and different 
fixed RBS dimensions {I = 5,10,15,20,25,30,35,40). Circles correspond to A = 0.5, dots - to A = 1 and 
crosses to A = 2. Red - method 2, blue - method 1. The y-axis is the decimal logarithm of the error. 
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Similarly, the error from Method 1 when using 20 snapshots (blue dots, Figure 6) is larger or equal to 
the error from Method 2 (red dots) when the dimension of the RBS is greater than 20 (third row of plots). 
In Figure S6 (Supplement) we present the distributions of the singular values in the 6 calculations (Method 
1 and Method 2, three valued of A). 

These results demonstrate the possibility of reducing the error of approximation using ROMs constructed 
via Method 2 compared to ROMs via Method 1 while keeping a low dimension of the ROM. 

5 Discussion 

Time derivatives include important information about the behavior of solutions of time dependent problems. 
It is reasonable to expect that POD ROMs including this additional information would approximate more 
accurately the solution of the FOM. It can be argued that difference quotients, used by some authors as 
discussed in the Introduction, are good approximations of the time-derivatives, which is true when the time 
intervals between snapshots are sufficiently small. DQs are not approximations if the solution snapshots are 
taken at relatively large intervals. Further, for dynamical systems of the form (HID calculating the derivatives 
at the times at which solution snapshots were calculated, come at almost no additional computational cost. 
Finally, unlike DQs, derivative snapshots do not belong to the reduced space generated by the solution 
snapshots. These arguments justify a comparative study of the accuracy of the POD ROM based on two 
types of snapshot selection - with and without using time derivatives as snapshots. 

We have derived error bounds that suggest that using time derivative snapshots may decrease the ap¬ 
proximation error if the first neglected singular value is not too large. We also demonstrate by numerical 
examples that the method with time derivative snapshots can yield significantly smaller (error of approxi¬ 
mation. Specifically, we show that when we take a small number of snapshots and use all of them to define 
the ROMs, Method 2 produces considerably smaller approximation error than Method 1. 

The general method we use to derive the error bounds in this paper has been used by other authors: 
deriving an equation for the error, directly integrating it and applying the Gronwall inequality by using 
assumptions for sufficient smoothness of the right-hand side of the system and its solutions and the respective 
Lipschitz constants. The innovative part of the method used here is the application of interpolation methods 
(Lagrange and Hermite interpolations in the two cases considered). This approach enabled the derivation 
of error bounds containing the time intervals between the snapshots and the orders of approximation 
expressed in terms of A^. 

This analysis of the error is the first emphasizing the relative significance of the two sources of the error 
- one coming from the reduced dimension Z, via the size of the first neglected singular value of the snapshot 
matrix, and the other from the term 0(A“), where a has different values (4 and 2) for the two methods 
(with and without time derivative snapshots). We believe that the study presented here contributes to 
understanding the error in POD ROMs. It is notable that even though we derive upper bounds for the 
error, and not estimates, these bounds are quite informative about these sources of error and the numerical 
experiments support what is expected from the bounds. 

Based on the numerical experiments performed it seems that the order a of the error may be higher than 
predicted by the error bounds we derive. More accurate bounds result from using the logarithmic norm as 
in |18] (which would primarily affect the exponential term in the error bounds) or other methods of error 
estimation that yet need to be defined. Our further goal is to devise methods using these or improved bounds 
for rational selection of the time moments of the snapshots. 
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7 Supplement 


V{x,t): MOL with n= 402 eqs 
BCy0t=-r0:y1t=-r1, r0=1.00 r1=5.00, 
a= 0.00 D1= 15.00 D2= 10.00 lambda= 0.00 gamma= 5.00 mu= 10.00 



Figure SI: Experiment A. Solution of the full system, x-axis - time, y-axis - space. The plot represents a 
set of 201 plots of the calculated solution Vi{t), i=0, 200. 


Singular values: method 1 (.) and method 2 (—) 

BC y0t=-r0; y1t=-r1, r0=1.00 r1=5.00, x= 10.000 
a= 0.00 D1= 15.00 D2= 10.00 lambda= 0.00 gamma= 5.00 mu= 10.00 



Figure S2: Experiment A. Distribution of the singular values for the two ROMs and three different snapshot 
spacings. Dots: Method 1; continuous line: Method 2; Red: A = 0.0025; Blue: A = 0.005; Green: A = 0.01. 


22 






V{x,t): MOL with n= 402 eqs 

Periodic BC y0t=-r0*(sin(t)f, y1t=-r1*{sin(t))^, r0= 1.50 r1= 0.50, 
a=0.10 D1=5.00 D2=1 .o 6 lambda= 2.00 gamma=5.00 mu=1.00 



Figure S3: Experiment B. Solution of the full system, x-axis - time, y-axis - space. The plot represents a 
set of 201 plots of the calculated solution Vj{t), i=0, 200. 


Singular values: method 1 (.) and method 2 (—) 

Periodic BC y0t=-r0'(sin(t)f, y1t=-r1*{sin(t))^, r0= 1.50 r1= 0.50, x= 10.000 
a=0.10 D1=5.00 D2=1.00 lambda= 2.00 gamma= 5.00 mu=1.00 



Figure S4: Experiment B. Distribution of the singular values for the two ROMs and three different snapshot 
spacings. Dots: Method 1; continuous line: Method 2; Red: A = 0.01; Blue: A = 0.02; Green: A = 0.04. 
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V(x,t): MOL with n= 402 eqs 

Periodic BC y0t=-r0*(sin(t)f, y1t=-r1*(sin(t))^ r0=1.50 r1=0.50, 
a=0.10 D1=5.00 D2=1.00 lambda=2.00 gamma= 5.00 mij=1.00 



Figure S5: Experiment C. Solution of the full system, x-axis - time, y-axis - space. The plot represents a 
set of 201 plots of the calculated solution Vi{t), i=0, 200. 


Singular values: method 1 (.) and method 2 (—) 

Periodic BC y0t=-r0*(sin(t))^, y1 t=-r1 *(sin(t))^, r0= 1.50 r1 = 0.50, x= 10.000 
a=0.10 D1=5.00 D2= 1.00 lambda= 2.00 gamma= 5.00 mu=1.00 



Figure S6: Experiment C. Distribution of the singular values for the two ROMs and three different snapshot 
spacings. Dots: Method 1; continuous line: Method 2; Red: A = 0.0025; Blue: A = 0.005; Green: A = 0.01. 
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